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Abstract 

We use a Green's function method with Random Phase Approximation to show how magnetic 
correlations may affect electric polarization in multiferroic materials with magnetic-exchange-type 
magnetoelectric coupling. We use a model spin ^ ferromagnetic ferroelectric system but our re- 
sults are expected to apply to multiferroic materials with more complex magnetic structures. In 
particular, we find that transverse magnetic correlations result in a change in the free energy of 
the ferroelectric solutions leading to the possibility for thermal hysteresis of the electric polariza- 
tion above the magnetic Curie temperature. Although we are motivated by multiferroic materials, 
this problem represents a more general calculation of the effect of fluctuations on coupled order 
parameters. 

PACS numbers: 75.80. +q, 77.80.Dj, 75.30.Et 
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I. INTRODUCTION 



Early theories that calculated the effect of magnetoelectric (ME) coupling on the sponta- 
neous magnetization M and electricpolarization P in so-called multiferroic materials used 
a Landau free energy formalism l|, |2j, y|. All terms in the free energy are written as prod- 
ucts of the order parameters M and P and must satisfy the symmetry of the paraphase. 
For example, for a centrosymmetric paraphase, the lowest order magnetoelectric free en- 
ergy term is P 2 M 2 . Recently, more exotic ME coupling terms have been proposed which 
also represent free energy invariants, such as a Dzyaloshinskii-Moriya coupling in an anti- 
ferromagnet P ■ (M x L) 0, S], and a spin density wave coupling proposed by Betouras 
P ■ (7 V (M 2 ) + 7' [M (V ■ M) — (M ■ V) M] + . . .) [6]. These coupling terms explain the 
multiferroicity of new candidate materials. 

However, by writing the ME coupling in terms of the order parameters, the important 
contribution to the coupling by correlations is ignored. Most obviously, above the lowest 
ordering temperature (out of the magnetic and the ferroelectric Curie temperatures, T C M 
and T c p respectively) there can be no ME coupling effects predicted by a Landau theory in 
zero applied field. In magneto-dielectric materials, magnetic correlations have been shown 
to have an important effect on dielectric constants Q] . In this paper we theoretically study 
the effect that magnetic correlations may have on spontaneous electric polarization through 
ME coupling. We find that the electric ordering temperature is shifted, even when it is 
above the ferromagnetic ordering temperature. We also discover the possibility for thermal 
hysteresis in P for multiferroic systems when a particular ME coupling term is allowec . 
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While calculations published in the last few years go beyond Landau theory 0, Q, 111 
they fail to identify this effect. 

The more general problem of understanding how correlations affect second order phase 
transitions in systems with coupled order parameters has been approached in the past jl3], 
particularly using techniques from renormalization group theory 14j. For a one dimensional 
system, the coupled parameters can be calculated exactly [15[ but for more general systems 
this is not the case. Our Green's function technique, which in this paper only includes trans- 
verse magnetic correlations, represents a unique and new approach. It also gives meaningful 
results for all temperatures, unlike the renormalization group method which can only give 
information on the critical behavior. 
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In Sec. [TT] we detail the Green's function technique with Random Phase Approximation 
and derive a free energy for a ferromagnetic ferroelectric bulk material, which can be utilized 
to solve for the coupled order parameters M and P, as well as other thermodynamic quanti- 
ties. In Sec. IIHI we provide some results for ME coupling which is both linear and quadratic 
in the order parameter P. We show how the coupling linear in P leads to the possibility for 
thermal hysteresis by altering the free energy of different local energy minima. In Sec. II VI we 
summarize the results, discuss the limitations of the current theory and provide an outlook 
on future work. 



II. GREEN'S FUNCTION METHOD 



A model ferromagnetic ferroelectric system with so-called "isotropic" or "exchange" ME 
coupling [3j is treated. The methods presented can be extended to treat multi-sublattice 
magnets but become much more complicated. We assume that TL < T? since this is the 

n 

case for the majority of multiferroic materials [3fl and also since then the effects of magnetic 
correlations on P will be more significant. 

The ferroelectric system is modeled using Landau theory for second-order phase transi- 
tions with free energy density given by: 

F FE = l -aP 2 + - A PP* - EP, (1) 

where E is an applied electric field parallel to the spontaneous polarization, a = Ak B (T—T c ), 
A and (3 are phenomenological constants, ks is Boltzmann's constant and T is temperature. 
T c is the ferroelectric Curie temperature in the absence of ME coupling. Landau theory is 
valid near phase transitions so we assume that T C M and T c p are sufficiently close for the results 
to be valid. P is treated as a scalar quantity. Its direction relative to the magnetization is 
not relevant in this simple model. 

We aim to write a free energy for the magnetic system with ME coupling to combine 
with Eq. (OQ) in order to solve for both order parameters simultaneously. We start with a 
microscopic Hamiltonian and use a Green's function technique (GFT) with a Random Phase 
Approximation (RPA) to derive the free energy. 
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The starting spin Hamiltonian is 

H = -/ i ^^-I(j + rP + 7 P 2 )^5,-5 J 
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+ S*(S Z ) + (S z ) 2 ) , (2) 

where the first term is the Zeeman interaction with h = g^sHo and Hq is an applied 
magnetic field. The second term represents the exchange interaction with a Taylor series 
expansion for the weak contribution from electric polarization P. Sushkov et al. recently 
used a Hamiltonian of similar form to describe how the magnetic exchange interactions 
together with magnetostriction in i?Mn 2 Os (i?=Y,Bi) couple strongly the mag netic system 



to a soft phonon mode associated with a spontaneous electric polarization [16j. The sum is 
over nearest neighbor pairs of spins at sites i and j and S i = Sf ±iSf . The constants T and 
7 describe the ME coupling strength that is linear and quadratic respectively in the spatially 
averaged order parameter P. We assume for a ferromagnet that all sites i are equivalent. 

Longitudinal spin terms SfS* are ignored in Eq. (j2J), but transverse terms SfS~ are 
kept. The reason for ignoring the l ong itudinal correlations is that RPA is known to produce 
spurious solutions for (S*S*) [l7, 18]. Although this does not cause significant problems 



when calculating the magnetization M = (S z ), it does introduce significant error in the 
subsequent calculation of the free energy. More advanced decoupling procedures may be used 
to gain better results for the longitudinal correlations [l^] but here we demonsrate possible 
effects of including transverse correlations as a first step. It should be noted that Callen 
decoupling [l9j] has the same problems as RPA with regard to longitudinal correlations and 
was specifically designed to improve RPA only when treating weak single-site anisotropies. 
We ignore anisotropies here since we will do the calculation for a spin | system. 



We define retarded Green's functions 
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G^t) = ((S+;S-)) = -i0(t)([St(t),Sr}), (3) 

where the square brackets indicate a commutator such that [A, B] = AB — BA, the single 
angled brackets indicate a statistical thermal average and the 9(t) function is a unit step 
function. The time Fourier transform of Eq. ([3]) is given by 

1 f +OD 

G^u) = ((St; 5J-» W = — / dt G«(t)e-^. (4) 

./ _ 



Then the equation of motion for the time Fourier-transformed Green's function is 

wC?«(«) = 6 f([S\S-]} + (([Sf,H];S.)U (5) 
where 5ij is the discrete Kronecker delta function. Substituting Eq. (j2J) into Eq. (jSJ) gives 

= ^(25*) (6) 



+ (/i + ^(j + rp + 7 p 2 )(Sf))G 



-(J + rP + 7 P 2 )^((^ + ;57)) w , 
z 

where 2; is the number of nearest neighbors to a site. The sum over I is over the nearest 
neighbors to site i. 

The last term in Eq. ([6]) represents a higher order Green's function which can be approx- 
imated using RPA: 

((S?S+; Sr))„ ~ (S*)«5+; 3T» U = (£*)Gy(w) (7) 
in order to obtain a solution for Gij(u). We perform a spatial Fourier transform 

G(w,fc) = G l3 {uoy~ lHri ~ rj) (8) 

= ^Y, G ^ k y Hr '~ rj) > ( 9 ) 

where the sum over 77 — 77 is over all displacements from site i to site j, and solve Eq. (jBJ): 



G(u,k) = . {SZ) (10) 
7r(a; - to k ) 

uj k = (J + TP + 1 P 2 )z{S z )(l- lk ) + h. (11) 

The structure factor 7^ = - £\ e * fe -°» ( where otj is the displacement from the reference site 
to its neighbor i, should not be confused with the coefficient for ME coupling 7. 

Eq. (TTTT) gives the dispersion relation for magnons within the RPA. The average electric 
polarization can be seen to alter the dispersion and application of an electric field will 
be able to shift the frequencies via the isotropic ME coupling. The possibility to tune 
spin wave dispersion using electric fields may have potential application in spin wave logic 



devices 
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231 ] . If our model allowed for fluctuations in the electric polarization, rather 
than just the magnetization, then the resonant modes may be "electromagnons" with a dual 
magnetic/electric nature. 



The transverse correlation function between neighboring spins i and j is calculated from 



Eq. (TTU1) using the Spectral Theorem 



21j with the aid of complex variable methods: 



t^O k 



(sjsf) = ^^J2 e ~ lk ' [ri ~ rj) / ^— — 



x 



7r(u; + — cjfcj 7r(u; — «e — Wfc) 
1 >^ 2(S z )e- ik < ri ~ r ^ 



fe_ 



(12) 



For spin | the magnetization is given by |2y, ]24| 



(13) 



M= (S z ) = i(l + 2$)- 1 , 
where $ = Sfc( efcsT ~~ -Q _1 - For general spin the result is 

M _ (s - m + ®r +1 + (s + 1 + m 2S+i (u] 

(1 + $)25+l _ $25+1 • ^ 

Having found a self-consistent equation for M, we need also to find an equation for P 



25| to derive the free 



using the free energy. We follow the workings in Appendix A of Ref. 
energy. 

The expectation value of the single-site Hamiltonian, derived from Eqs. ([2]) and ( TT21) . 
gives the intrinsic energy $ per magnetic lattice site: 



£ = (Hi) 



Ms z )--(J + tp + 1 p 2 )(s z ) 



x 



Tfe 



N V(e^ 



(15) 



From the intrinsic energy, we can derive an expression for the free energy F by making use 
of the relations: 



F = <S-TS, 
S= 

\dT /M 

where 5* is the entropy. Rearranging these we obtain 

F(r) = <f (o)-r/ r d / (r) -/ (0) . 



(16) 
(17) 

(18) 



This free energy is not a function of M = (S z ) since M must be constant in the definition 
of S which we use [Eq. ( 117]) ]. Substituting Eq. ([15]) into Eq. (ITS]) we obtain the free energy. 



Adding it to the ferroelectric free energy [Eq. (pQ)], the total free energy per magnetic unit 
cell of volume V is given by 



GFT 



v(\aI*+\pP i -Ep\ (19) 
-h(S z ) - |(J + TP + 1 P 2 )(S Z ) 

N ^ U k \ 
k 



X 



This expression is true for general spin. The last term in Eq. ([19]) represents the contribution 
from the transverse correlations and at low temperatures gives the free energy of a gas of 
noninteracting magnons, proportional to fc^T^ fc ln(l — e k B T ). A solution for P can be 
found by numerically minimizing Eq. (TT9~D . 

To find M(T) and P{T) we use an iterative procedure, starting at low temperatures 
and using Mq = ±| as an initial point. Substituting Mo into Eq. ffT9l) and minimizing, we 
obtain P x . Substituting P\ into the right hand side of Eq. (|T3|) . we obtain M\. This process 
is repeated n times until there is no longer a change in M n and P n to the required precision. 
For higher temperatures, the most useful starting point for iteration is the solution already 
found for lower temperatures. 

We will compare our results to those found using mean field theory (MFT) in order 
to examine the effect that the transverse magnetic correlations have on the solution. MFT 
corresponds to ignoring the transverse correlation term between neighbors, S^S~, in Eq. (T5]) 
and thus corresponds to reducing the problem to a single particle problem. The resulting 
free energy for spin | can be found simply using the magnetic partition function Z: 

Fmft = F FE - k B T In Z 

= F FE -k B T\nl (Si\e~^\Sf) 

= V QaP 2 + \(3P 4 - EP^ (20) 

- Z -(J + TP + 1 P 2 ){S Z ) 2 

h + z(S z }(J + TP + 7P 2 



— ksTln 



^cosh ^~ 



2k B T 

This free energy may be minimized with respect to both M and P, unlike Eq. (fT9|) . to solve 
for the order parameters at a given temperature T. 



III. RESULTS 



A. Changes to critical temperatures 

In Fig. [T] the magnetization (panel a) and electric polarization (panel b) are plotted as 
a function of normalized temperature IzbT/zJ for a, S = \ ferromagnetic ferroelectric with 
isotropic magnetoelectric coupling that is quadratic in P. The MFT results (Eq. (I20p ) 
are shown by the solid lines and the GFT results are shown by the dots. The material 
parameters used are a = Ak B (T — T c ), k B T c /zJ = 1/3, A = (3 = 1, T = and •y/zJ = 0.05. 
The coupling strength 7 is unphysically large but allows the effects of the isotropic coupling 
to be seen clearly. Also, a small applied field given by h/zJ = 0.0017 is applied in order 
that the magnetic Green's function is defined at all temperatures. 
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FIG. 1: (Color online) The magnetization (a) and electric polarization (b) plotted as a function of 
normalised temperature k^T j ' zJ of a ferromagnetic ferroelectric with isotropic ME coupling that is 
quadratic in P, calculated using mean field theory (solid lines) and using Green's function method 
(dots). The parameters used are a = Aks(T — T c ), ksT c /zJ = 1/3, A = (3 = 1, T = 0, 7 = 0.05zJ 
and h = 0.0017zJ. A simple cubic lattice is assumed. 

From Fig.dla), the GFT gives a value for the Curie temperature T C M lower than the mean 
field value. The reason is that an introduction of correlations results in less thermal energy 




S 



being necessary to flip spins, thereby lowering the critical temperature. From Fig. [T](b), 
we see that the Green's function with RPA treatment for the spin system gives a critical 
ferroelectric temperature which is also lower than the mean field prediction. Most stikingly, 
the magnetic transverse correlations affect the ferroelectric system above the magnetic Curie 
temperature and cause a reduction in the spontaneous polarization. 

In Fig.[T]only the solutions corresponding to M > and P > are shown. However, there 
are four solutions which correspond to local energy minima for the ferromagnetic ferroelectric 
system in the four different quadrants of the phase space given by {M, P}. Because of the 
symmetry breaking applied field h > 0, the two solutions with {M > 0, P > 0} and 

{M > 0, P < 0} (call them S ++ and S" H respectively) are degenerate and correspond to the 

lowest energy solutions for isotropic coupling that is quadratic in P. Applying a symmetry- 
breaking positive electric field breaks the degeneracy and makes S ++ the equilibrium solution 
for all temperatures. This is not the case when we consider a multiferroic system with 
isotropic magnetoelectric coupling that is linear in P. 

B. Thermal hysteresis 

In Fig. [2] we show two solutions S ++ (both M and P positive) and >S + _ (M positive 
and P negative) as a function of normalized temperature ksT/zJ when the ME coupling 
is linear in P. S ++ has magnetization plotted in panel (a) and polarization in panel (b). 

has magnetization plotted in panel (c) and polarization in panel (d). The solid lines 

show the MFT results and the dots show the GFT results. The parameters used are given 
in the figure caption. The coupling linear in P is extremely rare in ferromagnets since it 
relies on broken inversion symmetry. However, it may exist in frustrated spin structures, as 
are typical for multiferroic materials. 

The most interesting deviation from the MFT results is that S ++ no longer exists above a 
temperature of ksT/zJ ~ 0.25 when transverse correlations are included (see Fig. [2](a) and 
(b)). This local free energy minima ceases to exist and there is a discontinuous ferroelectric 
transition from one solution to the other. To explain this result, we need to consider the 
free energy of the solutions. 

In Fig. [3] we show the free energy of the two solutions illustrated in Fig. [2] as a function 
of normalized temperature. S ++ is shown by solid lines and S + - is shown by dashed lines. 
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k B T/zJ k B T/zJ 

FIG. 2: (Color online) Two solutions for {M, P} in a ferromagnetic ferroelectric with isotropic ME 
coupling that is linear in P, calculated using mean field theory (solid lines) and Green's function 
methods (dots). Panels (a) and (b) show the magnetization and polarization as a function of 
temperature corresponding to solution <S++, which has positive P. Panels (c) and (d) show the 

magnetization and polarization as a function of temperature corresponding to solution , which 

has negative P. The material parameters used are a = AksiT — T c ), kBT c /zJ = 1/3, A = f3 = 1, 
r = 0.05zJ, 7 = and h = 0.0017ZJ. A simple cubic lattice is assumed. 

The MFT free energies [Eq. ( 120]) ] of the two solutions are shown in red and converge at 
high temperatures. The GFT results [Eq. (TT9]) ] are shown in black. When correlations are 
ignored, the MFT results show that S ++ is always the lowest energy solution. This is because 
the coupling gives rise to an effective symmetry-breaking electric field E eS = Tz 2 y 2 > (see 
Eq. flZOP )- This is also the case in the GFT for low temperatures. However, at higher 
temperatures where the transverse correlations become large, the situation is different. 

Solutions S ++ and S + _ have different effective exchange interactions = J + TP. In 
consequence, the contributions to the free energy due to transverse correlations are different 
in each case. The energy associated with the correlations can be deduced by subtracting 
the MFT free energy from that found using the GFT (see Fig. [3]). Solution S + - has a lower 
exchange energy associated with correlations and also a lower Curie temperature than S ++ . 
This means that at ksT/zJ ~ 0.2 the lowest free energy state changes from S ++ to S^- 



10 




0.1 0.2 0.3 

k B T/ Z J 



FIG. 3: (Color online) The free energy corresponding to solutions (dashed lines) and S ++ 

(solid lines) for the ferromagnetic ferroelectric with isotropic ME coupling. The GFT predictions 
are shown in black and the MFT predictions are shown in red. Material parameters used are the 
same as those in Fig. [2l 

(see Fig. [3]). Also, for k B T/zJ > 0.25, the transverse correlations add an additional energy 
that makes S ++ unstable. 

An intriguing consequence is the possibility of thermal hysteresis. The state of the system 
depends not only on its temperature but also on the temperature it has had in the past. 
If the system is heated to above k B T / zJ ~ 0.25 then it must have solution S + _ since 
S ++ vanishes. If it is then cooled back down to below ksT/zJ ~ 0.2, the system may 
remain in solution S + _ even though it is now a metastable solution (see Fig. [3]). There 
is a probability that thermal fluctuations may take the system back to state S ++ with a 
characteristic time for switching r which depends on the temperature and the height of 
the energy barrier separating the solutions. If r is quite large then we have a long-lived 
two state system which may be implemented for information storage, provided the thermal 
hysteresis occurs in a region spanning room temperature. The application is to thermally 
assisted ferroelectric recording, where the coupling to magnetic correlations has enabled 
there to be thermal hysteresis of the electric polarisation P. As far as multiferroics are 
concerned, the only experimental evidence of thermal hysteresis in P that we are aware of 
as yet is at commensurate-incommensurate spin density wave transitions in DyMr^Os at 



low temperature 



26] 



We must question how our results would be changed if longitudinal correlations were 
included in the theory. We expect our result of thermal hysteresis to still hold because 
of the following argument. At very low temperatures, S ++ is the stable solution of the 
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system and both the MFT and GFT agree on this result (see Fig. [3]). At high temperatures 
T > T+ + ,T+_ and with h — > 0, (S?Sj) = (SfSp due to symmetry considerations. So 
we can expect that the difference in the free energy of the two solutions would have the 
same sign and would be larger if longitudinal correlations were included. Therefore, at 

high temperatures the result that is the stable solution seems robust. With the low 

temperature and high temperature results assumed correct, the system necessarily switches 
from one solution to another at intermediate temperatures. 



IV. CONCLUSION 



We have presented a model which goes beyond Landau theory of second order phase 
transitions to examine the thermodynamic properties of a system with coupled order pa- 
rameters. We use a Green's function theory with Random Phase Approximation to include 
the effect of transverse magnetic correlations on the free energy of a multiferroic system with 
magnetization M and polarization P. We find that the ferroelectric transition temperature 
is shifted as compared with mean field predictions that ignore magnetic correlations. This 
has implications when estimating ME coupling strength in experiments that measure a shift 
in ferroelectric polarization on application of a magnetic field at finite temperature. We also 
show that if a ME coupling that is linear in P is symmetry-allowed in a system, then there 
is the possibility for thermal hysteresis of P. Which free energy minimum represents the 
lowest energy state switches from one solution to another at finite temperature due to each 
solution having different exchange energy contributions. 

Our method has the advantage of being well-defined at all temperatures, unlike renormal- 
ization group methods. However, longitudinal magnetic correlations are ignored to simplify 
the theory and also any fluctuations in the ferroelectric system are ignored as a first approx- 
imation. More sophisticated decoupling procedures exist to include longitudinal magnetic 
correlations without creating spurious results for the free energy [171 ] . Also, it is possible to 
go beyond the Landau treatment of the ferroelectric system using, say, a Hubbard model 
which is suitable for displacive ferroelectrics (see for example Ref. 27l|). 

If the method were to be used to model a specific multiferroic material, then an extension 
to treat multi-sublattice magnets would be necessary. This would be possible using existing 
Green's function methods for antiferromagnets (see for example Ref. 28[ or 291] ). 
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